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Abstract —We present a new approach to simulating mixtures 
of gas and dust in smoothed particle hydrodynamics (SPH). We 
show how the two-fluid equations can be rewritten to describe 
a single-fluid ‘mixture’ moving with the barycentric velocity, 
with each particle carrying a dust fraction. We show how this 
formulation can be implemented in SPH while preserving the 
conservation properties (i.e. conservation of mass of each phase, 
momentum and energy). We also show that the method solves two 
key issues with the two fluid approach: it avoids over-damping 
of the mixture when the drag is strong and prevents a problem 
with dust particles becoming trapped below the resolution of the 
gas. 

We also show how the general one-fluid formulation can be 
simplifled in the limit of strong drag (i.e. small grains) to the 
usual SPH equations plus a diffusion equation for the evolution 
of the dust fraction that can be evolved explicitly and does 
not require any implicit timestepping. We present tests of the 
simplifled formulation showing that it is accurate in the small 
grain/strong drag limit. We discuss some of the issues we have 
had to solve while developing this method and Anally present a 
preliminary application to dust settling in protoplanetary discs. 


the particulate size is finite, giving rise to buoyancy forces. 
We neglect these since the volume occupied by dust grains 
is negligibly small in most astrophysical problems. We also 
neglect the thermal coupling between the phases. 

B. The stopping time 

The physics of drag is best quantified by the ‘stopping time’ 
t = _Md_ (6) 

which is the characteristic timescale on which the two fiuids 
dissipate their relative motion and stick together at a common 
velocity. It is inversely proportional to the strength of the 
mutual drag force. Short stopping times mean that the mixture 
moves together at a common velocity (the velocity of the 
bary centre), long stopping times mean that the two fiuids move 
separately, but intermediate stopping times produce the most 
dissipation, since drag acts to dissipate the relative motion. 


I. Introduction 

A. Gas/dust mixtures 


Multiphase flows where particulate matter is carried by 
a fiuid are crucial to many problems in astrophysics and 
engineering. Perhaps the simplest example is that of gas and 
dust, described by 
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where the subscripts g and d refer to the gas and dust, respec¬ 
tively, p and V refer to density and velocity and u and P are the 
specific thermal energy and pressure of the gas, respectively. 
The basic physics is that gas is affected by pressure where 
gas is not, and that the two fiuids are coupled by a drag term 
(with drag coefficient K derived from the microphysics of the 
interaction) which is dissipative. Additional terms arise when 


C. Gas/dust mixtures in SPH 

The traditional approach to gas/dust mixtures in Smoothed 
Particle Hydrodynamics (SPH) is to discretise these equations 
using a different set of particles for each phase of the mixture 
Q, giving ( O, Q) 
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where we use labels a,b and c to refer to gas particles; i, j and 
k to refer to dust particles and u is the number of dimensions. 
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Fig. 1. Linear waves in a dust-gas mixture (the ‘dustywave’ problem), 
showing the SPH two fluid solution with 2 x 100 particles (black circles=gas; 
open circles=dust) after 5.5 periods, compared to the analytic solution for the 
gas (solid red line) and dust (dashed red line). At low and intermediate drag 
the solution is accurate, but at strong drag the numerical solution because the 
short lengthscale separating the two fluids is not resolved. 


In recent work 0, 0 we made a number of key im¬ 
provements to the discretisation of this set of equation in 
SPH. We found a factor of 10 improvement in accuracy at 
no additional cost by employing a double-hump shaped kernel 
{Daj in equations [9] and [T0|), to compute the drag terms instead 
of the usual bell shaped kernel {Waj)- We also presented 
an improved implicit integration method and generalised the 
earlier methods of (h, a to spatially variable smoothing 
lengths. 

D. Two problems with two fluids 

1) Overdamping: The first problem we found with devel¬ 
oping algorithms for dust/gas mixtures was that there were 
few simple test problems which could be used to benchmark 
the numerical solution. This led one of us (GL) to derive the 
complete analytic solution for linear waves in the mixture, 
which we published in O. We found this immensely useful 
and enlightening, and indeed it revealed a rather fundamental 
limitation to the two fluid formulation. Figure [T] shows a 
typical SPH two-fiuid solution after 5 ^ wave periods, solving 
Equations ©-(O for a one dimensional linear-wave in an 
equal mixture of dust and gas while varying the drag parameter 
K, in each case compared to the analytic solution in both the 
gas (solid red line) and dust (dashed red line). 

The behaviour of the analytic solution is intuitive — when 
the drag is small the solution corresponds to an undamped 
sound wave propagating in the gas. At very strong drag the 
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Fig. 2. Gas (left) and dust (right) column density in a two-fluid simulation 
of material orbiting in a protoplanetary disc. Once the dust smoothing length 
becomes smaller than the typical gas smoothing length (solid red circle shows 
a 2for a representative gas particle) the dust particles become artiflcial 
‘trapped’ in high density rings, due to the lack of mutual repulsion between 
SPH dust particles. 


solution also corresponds to an undamped soundwave, with 
the only effect being a change to the effective sound speed 
because of the weight of the dust being carried along by the 
gas. Importantly it is only at intermediate drag (stopping times 
comparable to the wave period) that the solution should be 
strongly damped. 

The SPH solution, by contrast, shows a strongly damped 
solution at high drag. The reason is that to accurately capture 
the physics when the drag is strong, one must resolve the 
(very) short lengthscale separating the two fluids. Indeed, we 
show in ll2i| that using a very large number of particles (up 
to 10,000 in ID) does reproduce the analytic solution, but 
the resolution requirement is prohibitive (corresponding to 
h ^ tgCs, where Cg is the sound speed). The effect of under¬ 
resolving this length scale is to mimic the effect of a much 
larger separation length, giving a solution closer to that of 
an intermediate drag which is highly dissipative. This spatial 
resolution requirement is in addition to the usual stability 
constraint on the timestep At < tg from the drag terms. In 
other words, the two fluid method requires an inflnite number 
of particles and an inflnite number of timesteps to correctly 
resolve the limit of a perfectly coupled mixture. 

2) Artiflcial trapping of dust particles: Because they do 
not feel any mutual repulsion, dust particles can also become 
artificially ‘trapped’ in high density regions. We found this 
to occur whenever the dust collects on a scale smaller than 
the local smoothing length of gas particles. An example is 
shown in Figure |2l showing the dust particles in the centre 
of a protoplanetary disc simulation forming artificial ‘rings’ 
once the dust smoothing length is smaller than the typical 
gas resolution (shown by the red circle). This problem could 
be only partially mitigated by using the maximum smoothing 
length of the two fiuids in the drag terms (Equations [9] and 
M — the simulation shown in Figure [2] does this but we still 
found particle trapping to occur. 

Both of these issues motivated us to develop an alternative 
approach that correctly captures the limit of strong drag/short 
stopping time, which is the subject of the present contribution. 
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11. Dusty gas with one fluid 

From the DUSTYWAVE analytic solution one can intuitively 
understand the physics of the mixture in the limit of infinite 
drag/short stopping time. This is the limit in which the fiuid 
behaves as a single fluid mixture with a modified sound speed. 
This motivated us to rewrite the mathematical equations in 
such a way that the limit of tg 0 is obvious from the 
mathematics. 

A. One fluid to rule them all 

To do this, we simply employed a change of variables. 
Instead of solving for Pg, Vg and v^, in 0 we rewrote 
the equations in terms of the total density p, the dust fraction 
e = pd/p, the barycentric velocity v and the relative velocity 
between the fluids Av, where 


_ PgVg -h PdVd 

Pg + Pd 

(13) 

Av = Vd - Vg. 

(14) 


The original variables can be written in terms of the new 
variables according to 


Pg = (1 - e)p, (15) 

Pd = ep, (16) 

Vg = V - eAv, (17) 

Vd = V + (1 - e)Av. (18) 


The evolution equations in the new variables, written in 
Lagrangian form with the time derivative defined for a single 
fiuid moving with the barycentric velocity. 
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where f represents body forces such as a gravitational poten¬ 
tial. This change of variables entails no loss of information 
but reveals a great deal of the physics. 
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Fig. 3. Solution to the DUSTYWAVE problem using the general one fluid 
SPH method, showing gas and dust properties reconstructed from the SPH 
mixture particles (black points and circles) compared to the analytic solution 
(red solid and dashed lines for gas and dust, respectively). The SPH solution 
is in excellent agreement with the analytic solution in all cases, unlike for the 
two fluid method (Figure [fl. 


B. Physical interpretation 

First, we see that in the limit of Av ^ 0 and e ^ 0 
the equations reduce to the usual equations of fiuid dynamics, 
which we know can be solved using standard SPH techniques. 
Second, we see the physics of the drag clearly in Equation!^ 
— the pressure gradient acts as the ‘source term’ to grow 
the differential velocity, while the stopping time appears as 
the characteristic decay timescale. That is, in the absence of 
pressure gradients the relative velocity decays exponentially 
to zero on the stopping time. Furthermore terms such as the 
anisotropic pressure in Equation [22| that are second order in 
Av can be neglected when the differential velocity is small, as 
occurs when the drag is strong and the stopping time is short. 
So it is clear that a discretisation based on Equations [20]-[24| 
will reduce to a standard SPH formulation in the previously 
‘impossible’ limit of tg ^ 0- 

C. SPH discretisation of the general one fluid formulation 

In principle there is no loss of information in solving 
Equations l20l - [^ instead of [T]-[5] but the two approaches are 
conceptually very different. The key difference is that we 
must now consider only one set of particles moving with the 
barycentre of the fiuid, advected according to 

dx 


(25) 
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These particles represent neither gas nor dust, but rather the 
‘mixture’. The dust fraction e and likewise the differential 
velocity Av are now intrinsic properties of the mixture that 
are advected along with the particles. The evolution of u 
is modified to express the fact that the particles are no 
longer simply gas particles. One must also ensure that the 
conservation laws are satisfied — in particular that the mass 
of each phase is separately conserved, just as it would be with 
the two fluid approach. 

Nevertheless it is reasonably straightforward to write down 
the SPH discretisation in a way that does satisfy the con¬ 
servation properties. The resulting equations, ignoring for the 
moment artificial dissipation terms, are given by Q 
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Fig. 4. Solution to a dust-gas shock problem with strong drag, showing the 
solution with the one fluid method (black points) compared to the analytic 
solution in the limit where the stopping time is small (red line). 


now correctly reproduced in all cases, including the limits of 
both strong and weak drag (i.e. for both ts ^ oo and 4^0). 
It is intuitive why this should be the case, since the resolution 
problem with the two fluid method referred to in Section Il-Dll 
was related to the need to resolve the physical separation 
length between gas and dust particles. In the one fiuid method 
there is no physical ‘separation’ of the ‘gas’ particles from the 
‘dust’ particles — every particle carries information about both 
phases. (Very!) careful inspection of Figure [3] would reveal 
this, since to visualise the gas and dust ‘particles’ we have 
simply made two copies of the same set of mixture particles, 
one with the gas properties reconstructed using (fTSl) and ([TtI) 
and one with the dust properties reconstructed using (O and 

(HI. 

A further confirmation that the one fiuid method solves the 
resolution issue is shown in Figured showing the propagation 
of a shock in the mixture when the drag is strong. The solution 
in this case should be identical to a hydrodynamic shock but 
with the shock propagating at the modified sound speed due 
to the weight of the dust. It can be seen that the one fiuid 
method gives results in excellent agreement with the analytic 
solution (red line), whereas in la we found that around 10,000 
particles in one dimension were needed to produce the correct 
solution with the two fiuid method. 


The full set of equations including the modifications to artifi¬ 
cial viscosity necessary for shock capturing are given in Q. 

D. Tests of the general one fluid method 

Figure [3] shows the solution obtained on the dusty wave 
problem by solving Equations [26l - [30l instead of CHH One 
can see immediately that there is no overdamping of the fiuid 
when the drag is strong. Furthermore, the analytic solution is 


It is also intuitively obvious that the one fiuid approach also 
solves the dust trapping problem — there is no separate set 
of ‘dust’ particles, so they cannot become trapped below the 
resolution of the gas. Instead, the resolution of both phases 
is now tied to the density of the mixture rather than being 
separate for each phase. Thus by definition the gas and the dust 
components are resolved with the same resolution (smoothing) 
length. 
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E. Limitations of the general one fluid approach 

Despite the general nature of the one fluid formulation and 
the removal of the spatial resolution criterion at high drag, 
there remains the stability constraint on the timestep from the 
stopping time if explicit timestepping is used to solve (1^- 
(1^ . i.e. 

At < ts. (31) 

As discussed in implementing implicit timestepping to 
solve Equations and [30| when tg is small is considerably 
simpler than implementing an implicit method for the two 
fluid method (0,1141,1111). Nevertheless it represents an extra 
complication to the algorithm. 

A further limitation is that the single fluid description is 
a fluid description. In particular the velocity held is assumed 
to be single-valued. This assumption breaks down in the limit 
where there is little or no drag (tg ^ oo) and the dust particles 
act as a pressure-less independent set of particles, uncoupled 
from the gas. Eor astrophysics this occurs for large km-sized 
planetesimals in protoplanetary discs, which oscillate freely 
around the midplane and form structures which are supported 
by the velocity dispersion of their constituent particles. The 
velocity held in this case can be multi-valued which can 
be captured when representing the dust as a separate set of 
particles but not if represented as part of a mixture. 

This means that the one fluid formulation is mainly useful 
in the limit where the stopping time is short compared to 
other timescales in the simulation, which is the limit of 
small (cm sized and smaller in protoplanetary discs) grains 
in astrophysics. However, in this limit we can also derive a 
much simpler method. 


HI. The diffusion approximation for the mixture 


A. The terminal velocity approximation 


In the limit where tg is smaller than the computational 
timestep. Equations l20l - [^ can be simplifled by neglecting 
terms that are second order in Av 0. This is known as the 
Terminal velocity’ approximation (e.g. 0, nni), since the 
time-dependence in Av can be neglected and for the case of 
hydrodynamics we simply have 
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In this limit we simply recover the usual equations of hydro¬ 
dynamics supplemented by an evolution equation (1^ for the 
dust fraction and an additional term in the energy equation 


(|3^ . The pressure is also modifled by the presence of the 
dust since it depends on the gas density rather than the total 
density, giving for example 

P={‘y-l)pgU = {l-e){j-l)fm. (37) 

The above change to the pressure, and hence the sound speed, 
produces the zeroth order effect of a ‘heavy fluid’ seen in the 
analytic solution for the DUSTYWAVE (Eigs. [T] and[3l). 

If desired, one can simplify the formulation even further by 
making the approximation that 4 = 0. This leaves only the 
‘heavy fluid’ effect, giving precisely the limit of sound waves 
propagating at the modifled sound speed. This limit would also 
imply — from Equation [35] — a constant dust-to-gas ratio, 
which is what is commonly assumed when interpreting most 
astronomical observations of dust in space. 


B. SPH discretisation 

The SPH discretisation is also much simpler. Both ([33]) 
and ([34]) are identical to the usual hydrodynamic equations, 
so the only question is how to discretise ([35]) and the extra 
term in ([36]). There are two possible ways of discretising 
these terms; either make the terminal velocity approximation 
([32]) in ([^-([30b. or simply discretise the terms in ([35]) 
and ([36]) directly. We compare both approaches in m, but 
And only minor differences, and so here present only the 
‘direct second derivatives’ version since it is both simpler and 
computationally more efficient. 

Adopting a standard approach to writing second derivatives 
in SPH lfT2l-lfT4l we discretise dlS-® according to 
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where Dg = e^tg and the scalar kernel function Fgb is deflned 
in the usual manner such that VWg6 = Fab^ab and Fgb = 

l[Fab{hg) -h Fgb{hb)]. 

The second term in Equation [41] is one of the most bizarre 
SPH discretisations we have ever come across, but derives 
from the discretisation of goll combined with the need to 
conserve energy. A remarkable proof is given in the appendix 
of m that this is indeed a correct discretisation of the 
corresponding term in ([36]). 

As with the general one-fluid method, all of the conservation 
properties of the two fluid method are satisfled by dsn-dill). 
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Fig. 5. Test problem showing diffusion of the dust fraction in a three 
dimensional calculation. The plot shows the SPH solution of the diffusion 
equation 03 at various times, compared to the analytic solution(s) given by 
the solid red lines. 



Fig. 6. Convergence study on the dust diffusion problem, showing second 
order accuracy of the method in L 2 . 


That is, the mass of each phase as well as the total mass, 
momentum, angular momentum and energy are all conserved 
exactly. 


C Timestep constraints in the diffusion approximation 

In a medium with a uniform density and temperature (such 
that P = CgPg = Cg(l — e)p). Equation [35l becomes a simple 
diffusion equation for the dust fraction 

^ = V • (7?Ve), (42) 

where the diffusion coefficient is given by 

7? = etscl- (43) 


This implies a stability constraint of the form 


Af < 


etscj ■ 


(44) 


This constraint, while second order in the resolution length, 
presents the opposite constraint on the timestep compared to 
both the two fiuid method (Section U) and the general one 
fiuid method (Section HB. That is, the timestep is inversely 
proportional to the stopping time, whereas with the other 
approaches the timestep is proportional to tg. 

This implies that Equations dH-diiii require only explicit 
timestepping to correctly capture the limit of small grains/short 
stopping time/strong drag — a stark contrast to the situation 
with the two fiuid method. 


D. Tests of the diffusion method 

1) Diffusion of the dust fraction: The nature of the approx¬ 
imation suggests a simple test problem where the particles 
positions are fixed but the dust fraction is allowed to evolve 
according to (1^ . Assuming a pressure proportional to the gas 
density such that the evolution reduces to (|4^ . the analytic 


solution is well known, and hence can be used. Although 
the problem is not particularly physical (since the other 
fiuid properties are held fixed) it is useful to benchmark the 
discretisation of the diffusion equation. 

Eigure[5]shows the results of a three dimensional simulation 
employing 50^ particles where the dust fraction was initially 
set to 



with the numerical solution compared to the analytic solution 
at various times, assuming a constant tg = 0.1 and density 
and sound speed of unity. The evolution is indistinguishable 
from the analytic solution. This is quantified in Eigure [6] 
which shows the L 2 error between the numerical and analytic 
solutions as a function of resolution (particle spacing) for 
two different discretisations. As expected, the convergence 
is second order in the particle separation since the particles 
remain in an ordered configuration. 

2 ) Waves: The numerical results using our SPH implemen¬ 
tation of the diffusion method on the DUSTYWAVE problem are 
shown in Eigure |71 compared to the analytic solution given by 
the solid line. The damping in the mixture is correctly captured 
both when the drag is strong and even in the regime where the 
diffusion approximation is no longer valid. The small phase 
error in the solution with K = 10 is not from the method 
itself, but from an inconsistency in the initial conditions used 
to set up the problem (we start with Av = 0 but the diffusion 
approximation assumes that Av is finite at all times). 

3) Shocks: The results on the dusty shock problem obtained 
by solving (|38ll-(|4l]i are indistinguishable from that obtained 
with the general one fiuid method (Eigure |4]), with the advan¬ 
tage that only explicit timestepping is required. 

4) Fall of a layer of dust: A simple problem to illustrate 
that the diffusion method correctly captures the behaviour of 
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Fig. 7. SPH simulation of the dustywave problem using the diffusion 
approximation, compared to the analytic solutions given by the red solid and 
dashed lines (c.f. Figures \T\ and [3}. The diffusion approximation assumes 
that the stopping time is short compared to the timestep and is therefore 
only applicable when the drag is strong. However, it requires only explicit 
timestepping to correctly capture the limit of small grains. 


the dust is the ‘fall of a layer of dust’ problem introduced by 
Monaghan IS- The setup is a box with a vertical gravitational 
force, with the gas density stratified to balance the vertical 
gravity and periodic boundaries in the horizontal direction. 
We use a constant drag coefficient K = 10, an isothermal 
equation of state P = c^pg with Cg = 10 and a gas density 
Pg = 1 at ^ = 0. We set the density stratification, and in 
the one fiuid method the dust location, by slightly altering the 
particle mass. The dust is initially confined in 0.6 < p < 0.8 
with pd = 0.1 in the layer and 0 elsewhere. 


Figure [8] compares the results using the diffusion approx¬ 
imation (left) to the results obtained with the standard two 
fiuid method (right). The one fiuid method benefits from the 
regularisation of the mixture particles by the gas pressure, 
whereas discreteness effects are visible in the two fiuid method 
because the dust particles have no self-interaction. Neverthe¬ 
less, comparable solutions are obtained with both methods. 

5 ) Settling of dust in a protoplanetary disc: Our final test 
problem is drawn from our intended application, namely the 
settling and migration of dust in discs around young stars 
during the planet formation process. We consider a vertical 
section of disc at a particular radius, with the gas pressure 
in hydrostatic equilibrium with the vertical component of the 
gravity from the central star. We perform the test at 50 AU 
with mm-sized dust grains using an Epstein drag prescription 
(for full details see CD). The stopping time for grains of this 
size is a few percent of the orbital timescale, meaning that the 
terminal velocity approximation is valid. However the problem 
is still (just) tractable with the two fiuid method enabling us 
to compare all three approaches. 

Figure [9] compares the numerical results with the diffusion 
method (top) with the general one fiuid method (centre) and 
the two fiuid method (bottom). The left panel shows the gas 
density while the right panel shows the dust density after 20 
orbits. Settling of dust to the midplane is expected to occur 
on a timescale of ^ 10^ orbits. The dust density is better 
resolved with the two fiuid method because the resolution in 
the dust follows the dust mass (in contrast to being tied to 
the total mass in the one fiuid case) but it is also more noisy 
because there is no pressure force to regularise the dust particle 
distribution. Nevertheless the results demonstrate that the basic 
physics can be captured with any of the three approaches we 
have discussed in the paper. 

The solution with the diffusion method is ~50 times faster 
to compute, since it requires only explicit timestepping. 
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IV. Conclusion 

We have derived a completely new approach to simulating 
two-phase mixtures in SPH. Rather than modelling the mixture 
with two separate sets of particles, we have shown how the 
equations can be rewritten to describe a single fluid mixture. 
This solves two key problems with the two fluid approach in 
the limit where the drag stopping time is short compared to 
the computational timestep. Furthermore the equations in this 
limit reduce to simply the usual SPH equations supplemented 
by a diffusion equation for the dust fraction and a modiflcation 
to the energy equation and require only explicit timestepping 
in the limit of small grains. We have so far applied this to 
dust settling in protoplanetary discs but expect the method to 
be equally applicable for industrial and engineering problems 
where particulate matter is carried by a fluid. 
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Fig. 9. Settling of mm-sized dust grains in a vertical section of a protoplan¬ 
etary disc, showing gas density (left) and dust density (right) computed after 
20 orbits using the diffusion approximation (top), the full one fluid method 
(centre) and the two fluid method (bottom). The solution with the diffusion 
method is comparable to the other two approaches but is approximately 50 
times faster than the two fluid method to compute since it requires only explicit 
timestepping. 
























